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Computer simulations were used to study the gel transition occurring in colloidal systems with 
short range attractions. A colloid-polymer mixture was modelled and the results were compared with 
mode coupling theory expectations and with the results for other systems (hard spheres and Lennard 
Jones). The self-intermediate scattering function and the mean squared displacement were used as 
the main dynamical quantities. Two different colloid packing fractions have been studied. For 
the lower packing fraction, a-scaling holds and the wave- vector analysis of the correlation function 
shows that gelation is a regular non-ergodicity transition within MCT. The leading mechanism 
for this novel non-ergodicity transition is identified as bond formation caused by the short range 
attraction. The time scale and diffusion coefficient also show qualitatively the expected behaviour, 
although different exponents are found for the power-law divergences of these two quantities. The 
non-Gaussian parameter was also studied and very large correction to Gaussian behaviour found. 
The system with higher colloid packing fraction shows indications of a nearby high-order singularity, 
causing a-scaling to fail, but the general expectations for non-ergodicity transitions still hold. 



I. INTRODUCTION 

Colloidal suspensions are often referred to as model 
systems for studying fundamental problems in condensed 
matter physics Q . Most of the properties of colloidal sys- 
tems are similar to those of simple liquids, except for the 
difference in the time scales involved in the processes in 
liquids or colloids, making the latter more useful in the 
study of some basic questions. Moreover, the interaction 
forces between particles in a colloidal system are easily 
tailored (e.g. by adding salt or polymer). However, there 
are some features found only in colloids, such as aggrega- 
tion or gelation, which makes the study of these systems 
even more fascinating. 

Gel formation, or gelation, is found in systems with 
strong short-range attractions, and is a universal phe- 
nomenon observed experimentally in many different sys- 
tems, ranging from colloid-polymer mixtures |^| to 
charged systems 0, or to globular protein systems ||. 
Gelation is the formation of a percolating network (typ- 
ically fractal) of dense and more dilute regions of par- 
ticles with voids which coarsen up to a certain size and 
freeze when the gel is formed. This process is observed 
in the structure factor as a low-q scattering peak which 
moves to lower q, increasing its height, and then arrests 
H [?], ||. Description of this phenomenon has been at- 
tempted with percolation theories, theories of phase sep- 
aration for states inside the liquid-gas binodal (which 
is meta-stable with respect to fluid-solid coexistence for 
short interaction ranges) or in terms of a glass transition 
of cluster of particles [|[ |) . 

Recently, acknowledging its non-equilibrium charac- 



ter, gelation has been interpreted using the formalism 
of mode coupling theory, MCT, for non-ergodicity tran- 
sitions [h], |ll[ Q. This approach views the gel as par- 
ticles trapped by a network of bonds which hinders the 
particle motion, resulting in a non-ergodic state. Thus, 
gelation is caused by formation of long lived bonds, whose 
collective arrest is described clS cL cL normal non-ergodicity 
transition. (This is distinct from many earlier approaches 
whereby the bonds were assumed to form irreversibly 
from the outset.) In the present simulation study, we 
want to test this suggestion critically, thereby establish- 
ing the existence or otherwise of a non-ergodicity transi- 
tion corresponding to bonding network formation. 

Also present in colloidal systems is the equivalent of the 
usual glass transition in simple liquids, which occurs at 
high densities, and is driven by steric imprisonment. This 
transition has been studied experimentally and compared 
to MCT thoroughly |||, (h| g|, f(|. When two different 
non-ergodicity transitions are observed in a system, MCT 
predicts a high order singularity in the region where the 
driving mechanisms for both transitions are present |l7], 
|l8| , |l9f| . Therefore, a higher order transition is expected 
at high attraction strength and high density in colloidal 
systems with attractive interactions . 

Computer simulations have been used to test the ex- 
pectations from MCT in many different systems, such as 
a Lennard- Jones liquid |p2|, |25|, p4| p5[ , water |26|, [27|, 
strong glass formers |3^, p3[ p4| and polymers ||28| , ]29| , 
|30| , |3l| ] . The tests have shown that the predictions from 
MCT are correct, not only qualitatively but also, in part, 
quantitatively |33|, |35|. However, they have also pointed 
out some differences, especially in the spatial correlations 
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of particle mobility ||, ||, In none of these sim- 
ulated systems, however, did gelation occur, presumably 
because the attractions were not short-ranged enough. 

In this work, we have used molecular dynamics simu- 
lations to study the properties of the gel transition, and 
compared them with the predictions from MCT. (This 
was initiated in pc| where some further results may be 
found.) We take the numerous universal predictions of 
the theory to test the scenario qualitatively. Compar- 
ing with quantitative predictions available for systems of 
hard spheres 1 4l| ff2| ] , spheres with short range attractions 
]To| , |Tl] , p"2| , |2l]7~and the mentioned simulation studies, we 
identify the novel mechanism driving the non-ergodicity 
transition which is the cause of gelation for moderately 
dense suspensions. Molecular dynamics were used in- 
stead of Brownian dynamics because the choice of micro- 
scopic dynamics does not affect the relaxational dynam- 
ics of a system close to a non-ergodicity transition | p3 . 
By means of the Asakura-Oosawa interaction potential 
flf3l , we simulate the behaviour of a colloid-polymer mix- 



ture, which is a well- understood system 44 



For short interaction ranges, this system exhibits a fluid- 
crystal transition, at intermediate densities and increas- 
ing attraction strength, with a liquid-gas transition meta- 
stable to the fluid-crystal one. In our simulations, the 
system was modified to prevent both of these phase tran- 
sitions from occurring, in order to be able to study the 
transition from the fluid to the non-equilibrium states. 

The paper is organized as follows: Section II describes 
some results from MCT which will be used in the sub- 
sequent analysis of the simulation results. In section III 
the simulation method is presented and the details are 
given. Section IV, deals with the results and is divided 
into four subsections studying i) the correlation function, 
if) the time scale and the diffusion coefficient, Hi) the 
mean squared displacement and iv) a higher colloid con- 
centration. Finally, in section V, we present the conclu- 
sions of this work. 



II. MODE COUPLING THEORY 

In this section we will present the most important 
MCT results on non-ergodicity transitions. MCT at- 
tempts a description of the density correlator and its self 
part, in terms of a fluctuating- force correlator |4|, fl9[ . 
In this paper, only the self part of the density correlator 
will be studied, defined as: 



d>*(i) = (exptafoW-r^O))}) 



(1) 



where the brackets denote average over particle j and 
time origin, and q is the wave-vector. The equation of 
motion of <fr s q in Brownian (coarse grained) dynamics, is 
given by: 



Tgfit^t) + *£(<)+ / m q {t-t')d v %{t')dt' = 
Jo 



where r q is a single particle diffusive time scale and m q (t) 
is a mode coupling kernel which describes the cage ef- 
fect |3Cj. Within MCT, glass states are given by non- 
zero solutions of this equation for the long time limit of 
$>q(t — > oo) = f q , the so-called non-ergodicity parameter. 
It describes the glass structure and may also be called 
Lamb-Mdssbaucr factor. The glass transition is marked 
by a (generally) discontinuous transition from the unique 
trivial solution in the liquid, f q = 0, to multiple solutions 
in the glass, /| > 0, where only the highest solution is 
physical. Glass transitions can be classified according to 
the number, I — 1, of non-trivial solutions merging with 
the highest one, and the type of transition is noted as Ai . 

For liquid states close to the glass, a two step decay 
is observed for the correlator; the plateau is at f q and 
signals the proximity of the glass transition. Around this 
plateau, shows some universal properties, depending 
on the type of transition. For the most common type of 
transition, A2, the decay to the plateau, and that from 
the plateau, can both be expressed as power law expan- 
sions. In particular, the decay from the plateau is given 
by: 



*;(*) = f"-h q 1 Ht/r) b + h v 2 Ht/rf + 0((t/r) 3b ) (3) 



with hq 1 ^ and h q z> amplitudes and r the final or ct- 
relaxation time scale, b is known as the von Schwei- 
dler exponent, and depends on the details of the inter- 
action potential. Expression (||) implies time scaling for 
the decay from the plateau, called a-decay, for different 
states close to the glass transition. The time scale, r, 
diverges as the glass transition is approached according 
to a power law, with an exponent 7, which can be related 
to the von Schweidler exponent: r ~ |er| with a the 
distance to the transition |48, 49 1. On the other hand, 



,(2) 



the wave- vector dependence of the non-ergodicity param- 
eter and amplitudes gives some non-universal properties 
of the transition, providing information about the mech- 
anism causing the non-ergodicity transition. 

For high order singularities, the fluid states close by 
show again a two step decay in the correlation function, 
but the decays to and from the plateau are no longer 
power law expansions. Instead, logarithmic laws are ob- 
tained |l9). A salient feature is that a logarithmic 
decay around the plateau is predicted: 



= ft 



C,log(t/ti) 



(4) 



(2) 



where f q A is the non-ergodicity parameter of the high 
order singularity, C q is an amplitude and t\ is a time 
scale (the time when the correlator lies on the plateau). 

The mean squared displacement (MSD) can be stud- 
ied instead of the correlation function, obtaining a similar 
two step behaviour. Similar asymptotic laws to describe 
the decay to and from the plateau can be derived, and 
the parameters and exponents can be related to those of 
the correlation function W% . The value of the plateau in 
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the MSD defines the localization length and is a measure 
of the size of the cage. However, it should be noticed that 
the cage, as formed by other particles, is constantly re- 
structuring cooperatively. Only when the particles have 
broken free of their cages, diffusive motion is observed, 
with a self-diffusion coefficient, D s , that tends to zero as 
the glass transition is approached as D s ~ for the 
usual Ai transitions. 

Two different non-ergodicity transitions have been 
found in colloidal systems with a short range attraction 
]To| , : a steric hindrance driven glass transition and an 
attraction driven gel transition [llj [L2j. While the 
first is found at high densities and is qualitatively similar 
to the glass transition in the hard sphere system (HSS) or 
Lennard- Jones system (LJS), the gel transition occurs at 
high attraction strength for all volume fractions. Differ- 
ent properties for these two transitions are predicted, the 
main difference arising from the driving mechanism: the 
localization length is shorter in the gel than in the glass, 
resulting in higher non-ergodicity parameters. Also, a 
smaller von Schweidler exponent for the gel than for the 
glass is expected, implying a higher value of 7, i.e. the 
transition as observed by T q or D s is more abrupt. 

The actual shape of the non-ergodicity transition line 
depends on the details of the interaction potential, al- 
though some general features can be found. From lower 
to higher interaction strength, the glass line, is slanted 
to higher concentrations, showing that a weak attraction 
Quidizes the glass. However, at even higher interaction 
strengths, the gel transition occurs at lower colloid den- 
sity the higher the attraction strength. As a result, a 
re-entrance transition is obtained at high colloid volume 
fractions. The line may be wedge-shaped or curved in 
this region, depending on the range of the interaction. If 
the line is wedge-shaped a high-order transition (generi- 
cally A3) is present near the corner, whereas none exists 
if the line is smoothly continuous. An A4 singularity ap- 
pears right at the vanishing of the A3 point when the line 
first becomes smooth GIL 



III. SIMULATION DETAILS 

Equilibrium molecular dynamics simulations mimick- 
ing a colloid-polymer mixture were performed for a sys- 
tem composed of 1000 soft-core polydisperse colloidal 
particles. The core-core interaction between particles was 
modeled by: 



V sc {r) 



»12 



-36 



(5) 



where a\2 — a\ + 0,2, with a\ and are the radii of the 
interacting particles. A flat distribution of radii with a 
width of S = 0.1a, where a is the mean radius, was used. 
The exponent in V sc was selected high enough to avoid 
problems related to the softness of the potential |5ll] . The 



polymer induces an attractive depletion interaction be- 
tween the colloidal particles, which was modeled by the 
Asakura-Oosawa interaction potential [[l3[ |47|]. The ex- 
tension of this potential to take polydispersity into ac- 
count reads |52l: 



Vao{t) = -k B T(j)p 



fa " to) 



1 2 



(*?+!)- 



2i 



(6) 



for r < 2(ai2 + £) and for larger distances. Here, r\i = 
a i/£; V — fa + 772 ) / 2 , and 4> p is the volume fraction 
of the polymer. Note that the range of the potential is 
given by the polymer size, £, and its strength by <p p . This 
potential was modified around r — a±2 , to ensure that the 
minimum of the total potential (V sc + Vao) occurs at this 
point: for r < 2ai2+£/5 a parabolic form, which connects 
analytically to Vao at 2ai2 + £/5 and has a minimum 
in 2ai2, was used. In our simulations, the range of the 
interaction, 2£, was set to 0.2a, which would correspond 
to polymers with R g /a = 0.1 where R g is the radius of 
gyration. 

A long-range repulsive barrier was added to the inter- 
action potential in order to prevent liquid-gas separation 
(as shown below). The barrier had a maximal height of 
IksT, according to a fourth-order polynomial: 




for ro < r < r± and zero otherwise. The limits of the 
barrier were set to 7'o = 2(ai2 + £); an d r i = 4a, which 
was enough to prevent phase separation. The maximum 
height of the barrier equals the depth of the depletion in- 
teraction at contact for <j) p — 0.0625, much lower than the 
values where the gel transition takes place. The result- 
ing total interaction potential, Vtot — Vsc + Vko + Vbar, 
is analytical everywhere. It is shown in figure yj, where 
in order to indicate the spread induced by polydispersity, 
the potentials among three different pairs with differing 
radii are plotted. 

In our simulations, lengths were measured in units of 
the mean radius, a, and time in units of \J 4a 2 /3w 2 , where 
the thermal velocity v was set to -y/4/3. Equations of mo- 
tion were integrated using the velocity- Verlet algorithm, 
in the canonical ensemble (constant NTV), to mimick the 
colloidal dynamics. Every rit time steps, the velocity of 
the particles was re-scaled to assure constant tempera- 
ture. No effect of nt was observed for well equilibrated 
samples. The time step was set to 0.0025. Equilibration 
of the systems was tested by monitoring the total energy, 
and other order parameters (see below), and by measur- 
ing &g(t) and the MSD at different initial times. When 
the order parameters were constant and the &i(t) and 
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FIG. 1: Total pair interaction potential Vtot as function of 
the radial distance r = ri — ra| for three different particle 
pairs; a pair of particles with minimal radii a\ — <J2 = a — S, 
one with average radii a\ — ai = a, and one with maximal 
di = a 2 = a + 5 (from left to right). The inset shows the 
enlarged region of the attractive minimum. Crosses mark 
where the parabolic minimum smoothly matches to eq. (^|). 

MSD curves showed no dependency on the initial time 
(ageing), the system was considered to be equilibrated. 
The volume fraction of the colloidal particles, </> c = 

I™ 3 + (f ) ) n c , with n c the colloid number density, 
and the polymer volume fraction, cf> p , were the control 
parameters used to identify the states in the phase dia- 
gram. 

In order to explore the whole <f> p — <f> c plane in search of 
the gel transition, phase transitions which forbid access 
(in equilibrium) to important parts of the plane must be 
prevented. Several order parameters were used to iden- 
tify different kinds of ordering in our system and to mon- 
itor whether unwanted liquid-gas or fluid-crystal transi- 
tions were taking place. First, the onset of phase separa- 
tion involving states of different density can be detected 
by dividing the system into n 3 boxes and measuring the 
density in every box. The 'demixing' order parameter is 
defined as the standard deviation of the distribution of 
densities: 



*„ = £(p*-p) 2 (8) 

where pk is the density of particles in box k, and p is 
the mean density. This parameter is close to zero for 
an homogeneous system, and increases if it demixes into 
phases of different density. In our case, n has been set 
to 4, implying 64 boxes, and a box edge of about 5a 
(depending on <p c ). On the other hand, the orientational 
order parameter, Qq, as defined by Steinhardt et al. |5^, 
|54f , signals the presence of an ordered phase, and is used 
to detect crystallization. 

The phase diagram was probed using these parame- 
ters. In figure @ the results are presented for a bare 
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FIG. 2: Demixing (^4) and orientational (Qe) parameters for 
4> c — 0.40 and increasing polymer fraction, <j> p , for different 
systems: monodisperse without long-range barrier (squares), 
polydisperse without barrier (circles), and polydisperse with 
long-range barrier (crosses). 

system (monodisperse and without the long-range bar- 
rier), a polydisperse system without the long-range bar- 
rier, and the final system with both polydispersity and 
barrier. In this figure, the colloid volume fraction is con- 
stant, 4> c — 0.40, and the polymer concentration varies; 
an isochore is studied. The sudden increase in both ^4 
and Qq occurring at <p p — 0.20 for the bare system, sig- 
nals the crystallization boundary, in accordance with Di- 
jkstra et al. |l7J . Because of the short range of the poten- 
tial, this system has no liquid phase; i.e. the liquid-gas 
coexistence is meta-stable with respect to the crystal-gas 
transition. 

When polydispersity is introduced in the system, crys- 
tallization is prevented, as indicated by the constant 
trend of both parameters close to <p p = 0.20. However, as 
the system now does not crystallize, the liquid-gas tran- 
sition can be reached upon increasing the strength of the 
interaction. This demixing is signalled by an increase in 
^4, not involving local ordering. In order to avoid this 
separation, the long-range barrier has been introduced in 
the interaction potential. The energy of a dense phase is 
raised, and demixing is thus energetically unfavourable. 
Figure || shows that liquid-gas separation is indeed inhib- 
ited by the repulsive barrier. Instead, individual voids 
of finite size are created in the system, causing a low-q 
peak in the structure factor, S(q), presented in Figure 
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FIG. 3: Structure factors for different polymer fractions at 
4> c = 0.40: 4> p = gray line, 4> p = 0.2 dashed line, <j> p — 0.35 
dotted line and 4> p = 0.425 solid black line. Note the low-q 
peak raise as cf> p is increased. Inset: Pair distribution function, 
q(r) for the same states. 



|3|. There, S(q) is shown for different polymer fractions, 
ranging from no attraction {<p p = 0) to the closest state 
to the gel we have accessed {<j> p = 0.425). 

In the inset to Figure ||, the pair distribution function, 
g(r), is presented for the same states as the structure fac- 
tor. The value at contact, r = 2, increases continuously 
as the atraction strength grows, signalling increased lo- 
cal contact probabilities. This process will be shown to 
be responsible for the gel transition. In S q it becomes 
evident as an increase in the oscillations for large q. 

The low-q peak in the structure factor resembles the 
low-angle peak observed in light scattering experiments 
with colloidal gels || ||. However, whereas the peak 
in our system is an equilibrium property, induced by the 
specific shape of the interaction potential, the experimen- 
tal peak has a non-equilibrium origin. We also checked 
for the possibility of microphase separation, which in 
some cases can be induced by a repulsive barrier p5| . 
In our case, the small angle peak continuously increases 
with <pp, but stays finite and smaller than the neighbor- 
ing peak as we approach the gel transition. We interpret 
this to indicate that we do not have microphase separa- 
tion, and we also observed no other signs of such ordering. 
Furthermore, since the relevant wave vectors in the MCT 
calculation of the gel transition are the high ones (around 
2-7r/£), the change in low-q region in the structure factor 
is expected to have little effect on the gel transition. 



IV. RESULTS AND DISCUSSION 

This system has been previously shown to undergo 
both the glass and gel transitions as stated by MCT. 
It also exhibits a logarithmic decay in the correlation 
function at high colloid concentration, indicating a high- 
order singularity in that region flip] . In this section we 
will discuss the properties of the gel transition, and com- 



pare them with MCT and with those of the HSS and 
other systems, which are similar to the glass transition 
at high colloid concentration. We test for differences by 
comparing quantitatively the non-universal features of 
the transition, which will aid in the identification of the 
driving mechanism. 

The gel line is predicted to extend to low packing 
fractions with the same qualitative properties. In order 
to test these properties we have performed simulations 
at two different colloid concentrations, <fi c = 0.40 and 
(f) c = 0.50, where the gel line is far away from the per- 
colation one. At high concentration, the higher order 
singularity is expected to affect the equilibrium states, 
disturbing some features of the gel transition. 



A. Self-Intermediate Scattering Function 

The scaling prediction for the a-decay of states close 
to a non-ergodicity transition is tested in Figure ||, for 
constant colloid packing fraction, <fr c = 0.40. Two dif- 
ferent representative wave-vectors are presented in this 
Figure, q = 6.9 and q = 15. As observed at the glass 
transition in the HSS and many other different systems 

e Eai r 
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IfJ, the a-scaling property holds. In 
comparing these correlation functions with those typical 
for the HSS or LJS, it is noticed that in Figure [| the 
a-decay of the correlators is more stretched, implying 
a smaller von Schweidler exponent at the gel transition 
than at the glass transition. Because of this stretching in 
the a-decay, a clear plateau is not observed, although a 
slowing down of more than four decades is studied. Nev- 
ertheless, extrapolating the relaxation curves to extract 
plateau values, much higher values are found than at the 
glass transition in the HSS or the LJS. 

We have analyzed the state C = 0.40 and <j> p = 0.42 in 
more detail, which shows four decades of slowing down 
compared to the purely repulsive situation upon turning 
on the attraction. Because scaling is observed in Fig- 
ure [|, studying only one state is enough to analyse the 
a-decay of the correlation function. The slowest state, 
4> p = 0.425 was not chosen because it strongly deviates 
from the expected behaviour of r q vs. <p p (see Figure || 
and discussion thereafter). The correlation functions at 
different wave- vectors for state <j> c — 0.40 and (f> p = 0.42 
are presented in Figure |. The range of wave-vectors 
studied, where the plateau height changes is much wider 
than the range for a similar change in f q at the glass tran- 
sition of hard spheres or Lennard-Jones particles. This 
feature indicates that the relevant distances for the gel 
transition are much shorter than for the usual glass tran- 
sition. 

The correlation functions were measured until the av- 
erage particle displacement was 5a, which is one fourth of 
the box size (21.95 a). Thus, extending this measurement 
to longer times in order to observe the whole a-decay at 
low q is troublesome. If the diffusion coefficient diverges 
at the same rate as the a-time scale (as predicted by 
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FIG. 4: Intermediate scattering function (self part), vs. 
re-scaled time with the a-time scale, r q , for different states: 
<p c = 0.40 and P = 0.375, 0.39, 0.40, 0.41, 0.415, 0.42 and 
0.425 from right to left. Two different wave-vectors are stud- 
ied: q = 6.9 (upper panel) and q = 15 (lower panel), with the 
KWW fits (dashed line) included. 



MCT), this problem would not appear. Thus, we are 
also observing a discrepancy between both time scale di- 
vergences, that will be further discussed below. 

The impossibility of observing a clear plateau, as men- 
tioned above, makes it more difficult to analyse the cor- 
relators, since /* cannot be fixed a priori. Therefore, 
the cv-decay of the correlation functions has been fitted 
using expression (^|), with /*, h q and h q 2 ^ as fitting pa- 
rameters. The von Schweidler exponent was also fitted 
but was kept identical for different wave-vectors. It was 
found as b — 0.37, and the other results for the fitting 
parameters, are shown in Figure ^. The trends of these 
parameters are similar to that of the glass transitions in 
both HSS and LJS, but over a wider q-range in the gel 
case. This indicates that the localization length is quite 
different in the present system. The non-ergodicity pa- 
rameter exhibits a bell shaped curve, whereas the first 
order amplitude describes a maximum. The latter is de- 
termined from the fit up to a prefactor which depends on 
the choice of r in eq. (|3j). As an estimate we have used 
r q for q = 9.9 (r q is defined by ^ q {r q ) = f q /e), which 
yields values that are similar (in magnitude) to the HSS. 
The second amplitude shows a monotonously increasing 
behaviour with q, in accordance with the HSS, but it is 
always positive, unlike the HSS where it goes through 
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FIG. 5: Intermediate scattering function (self part), 
for the state C = 0.40 and cj> p = 0.42, 
for different wave-vectors; from top to botton: q — 
3.9, 6.9, 9.9, 15, 20. 25, 30, 35, 40, 50. The dashed lines are 
fittings from eq. (j3|) up to second order, with the same von 
Schweidler exponent for all q. 




FIG. 6: Non-ergodicity parameter, /| (open circles), and first 
(filled circles) and second order (diamonds) amplitudes from 
the fittings in Figure^] with b — 0.37 for all wave- vectors. The 
lines give the Gaussian approximation from eq. (^) fitted to 
all wave-vectors (solid line) and the three lowest wave- vectors 
(dashed line). 



zero at the peak of h q . 

The non-ergodicity parameter, can be approxi- 
mated using the Gaussian expression: 

f q « exp{-g 2 rf/6} (9) 

where r\ is the localization length. This approximation 
is known to be valid for low wave- vectors, and important 
deviations from the Gaussian behaviour are expected 
close to the glass transition. However, the value for the 
localization length obtained from fitting this curve (solid 
line in Fig. ^J), can be used as an estimate of the one in 
the MSD. 
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FIG. 7: Parameters used in the KWW fittings and compari- 
son with the von Schweidler fitting. Main figure: A q (closed 
circles) and /| (open circles), f3 q (crosses) and b (horizontal 
dashed line). Inset: T q (closed circles) and r q (open circles). 



The localization length so obtained is rf = 0.0126a 2 , 
much smaller than for the HSS or the LJS, where n is of 
the order of the Lindemann distance. This feature shows 
that the process causing the non-ergodic transition in our 
case has a typical distance much smaller than in the case 
of glass transitions in the HSS or the LJS. This agrees 
with the observation and discussion about the height of 
the plateaus, and of the different q-range covered by /* 
in Figure ^. Whereas the glass transition in the HSS is 
driven by core-core repulsions, the gel transition is caused 
by the short range attraction, therefore by bonds between 
particles, (see inset to Figure ||) whose size is of the order 
of the interaction range. An interesting analogy has been 
established between the mechanisms driving the forma- 
tion of gels and glasses, and the freezing transition |5(| . 

The a-decay of near-non-ergodic states can be also 
studied using the Kohlrausch- Williams- Watts (KWW) 
stretched exponential. The KWW expression is given 
by: 



Sf(t) = ^ g exp 



(10) 



where f3 q is known as Kohlrausch exponent, which has 
been shown to coincide with the von Schweidler expo- 
nent at high wave- vectors |57| . This expression has been 
fitted to very different systems, and describes the a-decay 
down to zero. We have fitted this expression to the a- 
decay in our system. However, since the correlators in 
Figure || do not show the complete a-decay, we have fit- 
ted expresion (10) to the master curve, obtained from the 
a-rescaling. Two of these fittings are presented in Figure 
^ by the dashed lines, showing that the KWW stretched 
exponential describes well the a-decay in this system. 

The fitting parameters, A qi (3 q and r g , are presented in 
Figure [?], and compared with the corresponding parame- 
ters in the von Schweidler formalism. In such a way, A q 



is compared with the non-ergodicity parameter, f3 q with 
the von Schweidler exponent, and the with r q . As 
expected, the height of the plateau can be determined 
equally well both by the KWW or von Schweidler anal- 
ysis. The same holds for the time scales, r q KWW and 
T q . The Kohlraush exponent is expected to tend to 1 
at low wave-vectors, and to approach the value of the 
von Schweidler exponent at high q. The low-g limit is 
explained because diffusion is the dominant process over 
long distances, whereas at short distances (comparable to 
the cage size) the dynamics is dominated by the cooper- 
ative local rearrangements. This behaviour is predicted 
from MCT J57j, and has been observed in different sys- 
tems, such as molecular glass formers pH , and in simu- 
lations of polymer melts |[l| , and of water . In our 
case, the low-g limit is not observed, but (3 rises as the 
wave-vector decreases, indicating that the expected be- 
havior may appear at lower q below the small angle peak 
in S(q). At high wave- vector, the Kohlraush exponent 
crosses the von Schweidler value, but stays close to it. 
Although an exact agreement is not observed, we may 
conclude that the correct general trend is obtained. 



B. Time Scale and Diffusion Coefficient 

An important universal prediction of MCT is the exis- 
tence of power law divergences for both the time scale, t, 
and the inverse of the self diffusion coefficient D s , with 
the same exponent in both cases, 7: 



,) 7 and D s ~ -4> p y (11) 



where (ftp is the polymer volume fraction where the gel 
transition occurs. The relation between exponent 7 and 
the von Schweidler exponent, 6, is also universally estab- 
lished by MCT g|. 

Testing of the power law divergence (and measuring 
of 7) is usually carried out plotting r q as a function of 
4>p — <ft p for different values of <frp, looking for a straight 
line. This method is cumbersome, even more as devia- 
tions from it are expected for states close to the transi- 
tion, and precise values for 7 and <ftp cannot be given. 
To avoid this difficulty, we have calculated 7 from 6, as 
given by MCT, and with this particular value of the ex- 
ponent looked for the power law divergence. In such a 
way, we are testing the compatibility of MCT predictions 
with our data. 

Figure || shows the wave-vector dependent time scale 
r q vs. 4>p — (ftp for different wave- vectors. For every wave- 
vector, T q ^ was extrapolated to zero, yielding a value for 
the polymer fraction at the gel transition, (ftp' q . The final 
value of (ftp, used in Figure ||, was calculated as the aver- 
age value for all wave-vectors studied. The linear trends 
in Figure || for (ftp — <fi p > 5 • 10~ 3 , shows the power- 
law behaviour predicted by MCT, with exponent 7 = 3.1 
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FIG. 8: Wave- vector dependent time scale, r q , vs. 4> p — 
(j> p for the isochore (j> c — 0.40 for different wave- 
vectors; symbols from top to botton correspond to q = 
3.9, 6.9, 9.9, 15, 20, 25, 30. The lines are power law fittings 
to q — 3.9 and q = 30. In all cases 7 keeps close to this two 
values. 



FIG. 9: Self diffusion coefficient, D s 



0.4519 (closed circles) and 



bp 



for two 
= 0.4265 



values of <f> p : 4 
(open circles). The dashed lines are the power law fittings to 
the data, with the exponent shown in the figure. Inset: D s T q 
vs. polymer fraction for different wave-vectors. Symbols as 
in Figure H. 



and tpp — 0.4265. The closest state to the gel transi- 
tion, 4>p — 0.425 deviates from the power-law behaviour 
observed for lower polymer fractions. Similar deviations 
have been observed in the HSS and LJS and can tenta- 
tively be attributed to thermally activated processes (or 
hopping events) M. 



As shown in eq. (11), MCT predicts a power-law for 
the self diffusion coefficient, D s , with the same exponent 
as the divergence of the time scale. Simulations on HSS 
and LJS have shown that a power law divergence is in- 
deed obtained, but with a different exponent than in the 
case of T q . Using the same procedure as described above 

(calculating 7 from 6 and extrapolating D s to obtain 



yields a value for 



0.4519, with the same 7 as 



for the time scale. This value of (j)^ is too far from that 
obtained using r q . Therefore, we cannot have similar (j)^ 
and 7 to explain the behaviour of both r q and D s , im- 



plying that the MCT prediction, eq. ( |11| ) is violated. 

In Figure we present D s vs. <fi p — 4> p using for <j>^ 
both the value estimated from r q and that from D s . We 
consider more desirable to have similar 0^ to explain the 
behaviour of r q and D s , even though this implies two 
different 7: 7 = 3.1 for r q and 7 = 1.23 for D s . As 
obtained in other non-ergodicity transitions |2^, [54| , 
the 7 exponent is lower in the diffusion coefficient than 
in the time scale, although the difference between both 
values of 7 is bigger in our case. 

In order to stress the different 7 exponents in the di- 
vergence of the time scale and 1/D S , we have plotted 
D s r q as a function of <j> p for different wave-vectors in the 
inset to figure ||. This product, that should be constant 
according to MCT, diverges as the polymer fraction ap- 
proaches 4>p ■ The divergence follows a power law with 
the exponent equal the difference between both values of 



7- 

The maximum in the self diffusion coefficient (upper- 
right corner of Figure ^) is a consequence of the re-entrant 
glass transition at high packing fractions [[n], ^| . A weak 
short range attraction at first destabilizes the cage and 
thus the glass transition moves to higher particle con- 
centration initially as the polymer fraction is increased. 
At constant colloid concentration the diffusion thus first 
speeds up with increasing <fi p , until for intermediate at- 
traction strengths the gel line is approached, where the 
opposite trend then dominates. At <fi c = 0.40, the glass 
transition is rather far removed and thus has little effect, 
but the increase of D is still measurable and the diffu- 
sion coefficient can be used as a measure of the distance 
to the closest transition. The maximum thus indicates 
the re-entrant shape of the non-ergodicity line. 

The wave-vector dependence of the time scale r q can 
also be compared with theoretical predictions. At low 
9, the time scale is expected to behave as q~ 2 , corre- 
sponding to a diffusive process over large distances. Yet, 
because the simulated scattering functions exhibit non- 
exponential relaxation even for the smallest wavevectors, 
this simple theoretical scenario is not expected to ap- 
pear in our case. At intermediate wave-vectors, where 
the Kohlrausch exponent becomes comparable to the von 
Schweidler one the theory predicts a decrease as q^ 1 ^. 
whereas at even higher q, the distances involved are dom- 
inated by the microscopic dynamics, and corrections to 
this behaviour are expected p3| . The inset to Figure [Toj 
shows T q for different states close to the gel transition. In 
order to make clear common properties the curves have 
been scaled vertically to collapse (main figure). 

It can be seen in this figure that the behaviour of r q at 
low wave- vectors (below q = 10), indeed shows a q~ 2 be- 
havior, which however is not the one explained by MCT. 
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FIG. 10: Inset: Time scale, r 9 as a function of the wave- 
vector, q, for different states close to the gel: <j> p = 0.42 (cir- 
cles), 4>p = 0.415 (squares), cj> p — 0.41 (diamonds), <f> p — 0.40 
(upward triangles), cj> p — 0.39 (left-ward triangles) and cj> p — 
0.375 (downward triangles) . Main body: Same data, re-scaled 
to collapse in the low-q power-law behavior. The dashed lines 
show power-law behaviours with exponents 2 (gray line) and 
3.3 (black line). 

At higher q, another power-law trend is observed, with a 
higher exponent: q~ 3 3 . The crossover from the low-g be- 
haviour to the high-q one, compares nicely with the wave- 
vector where the Kohlraush exponent becomes equal to 
the von Schweidler one (Figure [?]) . The exponent of the 
high-g region yields b = 0.30, lower than the value ob- 
tained from the analysis in Figure |5|. However, this value 
is quite close to the measured von Schweidler value and 
much smaller than the HSS one. Deviations from this 
power-law behaviour are observed at high q for the lowest 
4>p presented in the figure. These deviations are caused 
by the microscopic dynamics, as they occur when r q is 
lower than a certain value, regardless the polymer frac- 
tion. This value, presented in the inset as an horizontal 
line, is to ~ 0.6, which agrees with the time one would 
estimate from the correlators in Figure || 

C. Mean Squared Displacement 

We turn now our attention to the MSD curves, that 
were partially analysed to obtain the diffusion coefficients 
presented in Figure ||. We are only interested in the slow- 
ing down close to the gel transition and thus we do not 
show the MSD for low polymer fractions, where the at- 
traction speeds up the dynamics and increases the dif- 
fusivity (see figure 9 and H). The MSD, after a short 
initial regime of free flight, Sr 2 oc i 2 , slows down because 
of the particle interactions and takes longer and longer to 
reach the long-time regime diffusive, where Sr 2 — 6D s t. 
An important feature that can be obtained from the MSD 
of the particles in the system, is the localization length, 
where the particle interactions hinder particle motion 
most strongly, and in the idealized glass state, arrest it. 
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FIG. 11: Mean squared displacement of the particles as the 
gel transition is approached. Black curves from left to right: 
(p p = 0.30, 0.325, 0.35, 0.375, 0.39, 0.40, 0.41, 0.415, 0.42, 0.425. 
Grey curve: (f> p = 0. The dashed horizontal lines indicate 
the localization length corresponding to the glass transition 
(short dash) in the HSS and the estimate from Figure ^| (long 
dash). Dotted curve: mean squared displacement for a single 
particle in a frozen environment at (f> p = 0.425. 



It can be compared with the estimate using the Gaussian 
approximation (see Fig. ||). In Figure [Tl| we present the 
MSD for increasing polymer volume fractions. As the gel 
transition is approached, the localization length shows up 
as an indication of a plateau, signaling the bond forma- 
tion. As already discussed above, r/ is much shorter than 
in the HSS glass transition (upper dashed line in Figure 
|li"| ), because of the driving mechanism. 

The lower dashed line in this figure is the localization 
length, as estimated from the non-ergodicity parameter 
using the Gaussian approximation (rf = 0.0126). Al- 
though a clear plateau has not fully developed in our 
curves, its height seems to be above that estimate, by a 
factor ~ 1.5 — 2. Since the Gaussian approximation works 
very well in the case of the HSS, this suggests big non- 
Gaussian corrections at the gel transition. Before testing 
the Gaussian approximation, we stress that the local- 
ization length gives a typical size of the mesh of bonds 
formed between neighbouring particles, and that the slow 
structral units are continuously and cooperatively rear- 
ranging. In order to test this idea about a correlated 
region which cooperatively rearranges with and around 
each particle, a single mobile particle is considered in 
a fixed environment. A well equilibrated system with 
<f) p = 0.425 is frozen, and only one particle is allowed to 
move. This mobile particle now explores a frozen envi- 
ronment, providing the structural size of the region it is 
confined to. The mean squared displacement so obtained 
is given in Figure [ll] (dotted line). Some particles (1.6 %) 
were able to break their bonds and diffuse freely in the 
frozen environment. For the particles that stay localized, 
it can be observed that the length of the frozen bonds 
is much smaller than the localization length. This fact 
demonstrates that the structure of bonds, like the repul- 
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sive cage at the glass transition in the HSS or LJS, is 
dynamic, and constantly rearranges cooperatively. This 
collective restructuring of the system fluidizes it, and re- 
stores ergodicity, which cannot occur in the frozen sys- 
tem, where the particles are not able to diffuse even at 
very long times. 

We turn now back to the Gaussian approximation, and 
its accuracy. Usually, this is tested by measuring the non- 
Gaussian parameter, defined as: 



Q'2 



5(r 2 (i)) 2 



1 



(12) 



where the averages imply ensemble averaging. This pa- 
rameter measures the deviation of the probability density 
function for the single particle motion from Gaussian be- 
haviour, and vanishes for diffusive motion. Special care 
must be taken when performing the ensemble averages 
in polydisperse systems, as pointed out in p4j . The non- 
Gaussian parameter must be calculated for every particle 
(the averages in the definition above thus implying time- 
origin averaging only), and particle averaging is taken 
on the values of a.i (so long as long enough time inter- 
vals are studied, each particle will sample the distribution 
relevant to its own size in an ergodic fashion). The non- 
Gaussian parameters for states with increasing cf> p are 
presented as a function of time in Figure [l2|. At short 
times «2 tends to zero, since the system shows Gaus- 
sian behaviour during its unhindered ballistic regime. At 
long times, when the particles break free from their bonds 
and hydrodynamic diffusion holds, «2 again goes back to- 
wards zero. At intermediate times, corresponding to the 
plateaus in both the correlation function and the mean 
squared displacement, ct2 grows, since the single particle 
motion hindered by bonding is not Brownian. As a result, 
ai shows a maximum, whose height and position grows 
in time, because the particles take longer and longer to 
break free and start diffusing. 

The behaviour of the simulated ui obeys the general 
expectations |3(| |37], p9| , but important differences are 
observed in the comparison with the results for the HSS 
or LJS. Whereas in those cases the height of the max- 
imum for similar (or even higher) a-relaxation times is 
around 2, at the gel transition much higher values are 
measured. Another interesting difference is the failure 
of the short time scaling, observed both in the HSS and 
LJS. Both effects can be rationalized considering that 
the cage is indeed a network of bonds in the case of a 
gel, rather than a cavity. The strength of these bonds 
is given by the intensity of the interaction, and thus, it 
is modified for different states, disabling the short time 
collapse. Because the bonds are short ranged, they af- 
fect the particle motion from very short times onward, so 
that the particles feel the hindrance much longer in the 
gel case. 

It can be concluded that the non-Gaussian corrections 
are very important in the gel transition. Therefore, the 
localization length estimated from the non-ergodicity pa- 




FIG. 12: Non-Gaussian parameter 122 as a function of time 
for states approaching the gel transition at the same states as 
Figure [n]. The maximum increases with increasing (j> p . 



rameter may be inaccurate, as discussed above. However, 
it still provides an indication of how small the localization 
length is. A better indication of 77 can be obtained within 
the Gaussian approximation if only low wave- vectors are 
used in fitting expression (§. The fitted curve is pre- 
sented in Figure ^| by the dashed line, where only the 
three lowest q's are fitted. The estimated deviates 
from the data at higher wave- vectors, showing high non- 
Gaussian corrections. The localization length is higher 
than the previous value: rf = 0.0162. Thus, this fitting 
provides data more consistent with the MSD curves and 
the non-Gaussian parameter. 



D. Higher Colloid Volume Fraction 

We move now to a higher colloid volume fraction: 
4> c = 0.50. These results are presented to supplement 
the findings at the lower packing fraction and test for 
the prediction of stronger stretching closer to the higher 
order singularity. As indicated in the theoretical section, 
MCT predicts a higher order singularity in the vicinity 
of the junction of the gel and glass lines; i.e. at high 
polymer and colloid densities. In this particular system 
we found clear indications of this singularity in simula- 
tions at <j) c = 0.55 and cf> p — 0.375 pQ] . The isochore 
under study now, cj) c = 0.50, could be close enough to 
the higher order singularity to show some effects. 

In Figure [L5] we present the correlation functions for 
increasing polymer fractions at the same wave- vectors as 
Figure |[ re-scaled to collapse in the long-time decay. It is 
interesting to note that the polymer concentrations stud- 
ied in this case are lower than those studied at the lower 
colloid volume fraction. In accordance with experiments 
and theory, this indicates that the gel transition takes 
place at lower polymer fractions the higher the colloid 
concentration. 

In Figure 13, it can be observed that the correlators 
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FIG. 13: Correlation functions for <j> c — 0.50 and differ- 
ent polymer concentrations. ^From left to right: <j> p = 
0.35, 0.36, 0.375, 0.38, 0.385, 0.39. Two wave-vectors are 
studied, as labeled in the figures. The dashed lines represent 
the KWW fittings to the 4> c = 0.40 correlation functions. 

do not collapse over the whole a-decay, but only in the 
end. These deviations are expected because of the higher 
order singularity, which is at higher densities. However, 
we stress that although this singularity has clear effects 
on the correlation functions, they do not show so clear 
signatures as that of the <ft c — 0.55 isochore |l0|. At this 
high concentration, a logarithmic decay was observed, 
with a wave-vector dependent extension. 

Because these corrections affect the early a-decay, ana- 
lyzing the correlation functions is difficult. Furthermore, 
the plateau is not observed, and the von Schweidler anal- 
ysis is thus extremely difficult. In order to analyse the self 
intermediate scattering function, we compare the stretch- 
ing of the curves at <ft c = 0.40 and <ft c = 0.50; in Figure 
|l3| the a-decay master function of the <ft c — 0.40 state, 
as parametrized by the KWW fitting is included. It can 
be seen that this curve can be rescaled to collapse onto 
the a-decay of the correlators at <ft c = 0.50 for both wave- 
vectors at long times. This indicates that the von Schwei- 
dler exponent is very similar in both cases, but also points 
out the effect of the high order singularity. According to 
MCT, b should decrease as the singularity is approached, 
but this behaviour is not observed in our case. Com- 
parison of the a-decays by fitting the KWW stretched 
exponential to the master function is troublesome, since 
only the late decay is obtained unambiguosly. 

The similarity of both a-decays was used in the von 




io" 2 10'' 10° io 1 io 2 io 3 

t 



FIG. 14: Intermediate scatteringfunction (self part), <E>g, for 
the same wave- vectors as Figure |5| wave- vectors for <j> c — 0.50 
and (j> p — 0.39. The dashed lines are fittings from eq. Q up 
to second order, with the same von Schweidler exponent as 
4> c = 0.40. 



Schweidler analysis of the correlation function, and only 
the non-ergodicity parameter and amplitudes were fitted. 
Since the upper part of the decay is known to be affected 
by the higher order singularity close-by, that part must be 
discarded in the fittings. The correlation functions and 
fittings are presented in Figure |l~4| for the state <ft c = 0.50 
and (j) p = 0.39, for the same wave-vectors as Figure ||. 
The main conclusion is that the late a-decay at all wave- 
vectors can be correctly described by the von Schweidler 
decay, with the same exponent as the state at <ft c = 0.40. 
The non-ergodicity parameters obtained from the fitting 
are slightly lower than those of <ft c = 0.40, but similar 
within the error bars. According to MCT, decreases 
when approaching the glass part of the non-ergodicity 
line (signaling an increase in the localization length) . Our 
result is thus consistent with this prediction. 

With these values of the non-ergodicity parameter one 
can define also the wave- vector dependent time scale, r q , 
as discussed above. In order to test the value of the von 
Schweidler exponent, using eq. ( |ll|) we have performed a 
three parameter fitting to obtain 7 and ft^ . In Figure [l^, 
T q is presented as a function of — <ft p for different wave- 
vectors. The power-law fittings for two wave-vectors are 
also plotted and the critical polymer fraction (ftp is given. 

The values of 7 obtained from this analysis for differ- 
ent wave-vector range from 7 = 3.37 to 7 = 3.82, the 
mean value being 7 = 3.70. This value of 7 implies a 
smaller von Schweidler exponent, b = 0.33, in disagree- 
ment with our previous estimate, but backing the MCT 
prediction. Using the same value of (ftp , the vanishing 
of the self-diffusion coefficient, D s , can be analysed, and 
is presented in the inset to this figure. A power-law is 
observed in this case, with an exponent, 7 = 1.92, which, 
again in contradiction to MCT, leaves us with a big dif- 
ference between the two values of 7. 

The diffusion coefficients in the inset of figure Il5] again 
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FIG. 15: Wave- vector dependent time scale, r q , vs. <j)p—4> v for 
different wave-vectors; symbols from top to botton correspond 
to q = 3.9, 6.9, 9.9, 15, 20, 25, 30. The lines are power law 
fittings to q = 3.9 and q — 30. 4>? and 7 were fitted. 



indicate the re-entrat glass transition. They describe a 
maximum, more pronounced than that observed in Fig- 
ure H because the glass line is closer to the 4> c = 0.50 
isochore. The minimum in r q , which is observed only for 
q = 3.9, in an equivalent way indicates the shape of the 
non-ergodicity transition line. At higher wave-vectors, 
the glass transition causes very low / * and the time scales 
merge with the microscopic transient and thus this fea- 
ture is suppressed. 

The wave- vector dependence of r q can also be studied, 
as done for the lower concentration, yielding another es- 
timate of b. In this case, a similar plot as Figure [l^ is 
obtained, where the low q region is compatible with a q~ 2 
behaviour, and a higher exponent at higher q, yielding a 
value of b = 0.38. This value is in agreement with the 
nice comparison between the <j> c = 0.40 and 4> c = 0.50 
isochores, but not with 7 or the MCT prediction. We 
may then conclude that analysis of this state is extremely 
difficult, but our indications state that the von Schwei- 
dler exponent is similar for both packing fractions, but 
probably slightly lower in the higher concentration. 

Finally, we would like to point out that the non- 
Gaussian parameter at this packing fraction shows a be- 
haviour similar to that shown in Figure [l2], i.e. the peak 
is as high, and no short-time scaling is observed. 



V. CONCLUSIONS 

In this paper, by means of simulations, we have tested 
the universal predictions of MCT for gelation in colloidal 
systems, viewed as an attraction-driven glass transition. 
The self parts of the intermediate scattering function for 
states close to this transition have been analysed and the 



results were compared with the theoretical predictions. 
For the <fi c = 0.40 isochore, which is far enough from 
the high order singularity, the correlation functions can 
be cv-scaled. The time scale of the a-decay was shown 
to obey a power law divergence, with an exponent, 7, 
related to the von Schweidler exponent, obtained from 
the early a-decay. Both features are predicted by MCT 
for all non-ergodicity transitions. Also, the wave-vector 
analysis of the time scale follows the behaviours predicted 
by MCT, with a small difference in the value of the von 
Schweidler exponent. 

The wave-vector analysis of the correlation functions 
depends on details of the interaction potential, and thus 
provides information about the mechanism, leading to 
the transition. In our case, it establishes that the gel 
transition is driven by a short-range mechanism, namely, 
bond formation, as observed in the pair distribution func- 
tion. Additionally, it has been shown that the KWW 
stretched exponential can account for the a-decay of the 
correlation functions, as in other non-ergodicity transi- 
tions. 

We have also tested the Gaussian approximation, 
which works very well for the HSS. The non-Gaussian pa- 
rameter, «2, establishes that this approximation is much 
worse in the case of the gel transition than for the glass 
transition. It was also tested when comparing the es- 
timated localization length from the non-ergodicity pa- 
rameter with the MSD of the particles. The diffusion 
coefficient has been also studied. It tends to zero as the 
transition is approached following a power-law, with an 
exponent much lower than 7, in accordance with simu- 
lations of glass transitions in other systems, but in dis- 
agreement with MCT, where both exponents are equal. 

Finally, when the colloid concentration is increased, 
the system shows signatures of the high order singularity 
nearby and little can be discussed about the exponents 
b or 7. However, only slight changes in the numbers 
are expected, since the qualitative behaviour is repro- 
duced, except for the a-scaling. Also, the diffusion coef- 
ficient follows a power law with a different exponent and 
the non-Gaussian parameter reaches values similar to the 
cj) c = 0.40 case. 

Therefore, our main conclusion is that MCT accounts 
for most features of the simulated systems on approach 
to the gel transition, but the discrepancies already found 
in other non-ergodicity transitions (such as the repulsion- 
driven glass transition in hard sphere systems) are also 
obtained here. 
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